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Abstract 

We  consider  the  simulation  of  steady-state  variably  saturated  groundwater  flow  us¬ 
ing  Richards’  equation  (RE).  The  difficulties  associated  with  solving  RE  numerically 
are  well  known.  Most  discretization  approaches  for  RE  lead  to  nonlinear  systems 
that  are  large  and  difficult  to  solve.  The  solution  of  nonlinear  systems  for  steady- 
state  problems  can  be  particularly  challenging,  since  a  good  initial  guess  for  the 
steady-state  solution  is  often  hard  to  obtain,  and  the  resulting  linear  systems  may 
be  poorly  scaled.  Common  approaches  like  Picard  iteration  or  variations  of  New¬ 
ton’s  method  have  their  advantages  but  perform  poorly  with  standard  globalization 
techniques  under  certain  conditions. 

Pseudo-transient  continuation  has  been  used  in  computational  fluid  dynamics  for 
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some  time  to  obtain  steady-state  solutions  for  problems  in  which  Newton’s  method 
with  standard  line-search  strategies  fails.  It  combines  aspects  of  backward  Euler  time 
integration  and  Newton’s  method  to  select  intermediate  estimates  of  the  steady- 
state  solution.  Here,  we  examine  the  use  of  pseudo-transient  continuation  as  well  as 
Newton’s  method  combined  with  standard  globalization  techniques  for  steady-state 
problems  in  heterogeneous  domains.  We  investigate  the  methods’  performance  with 
direct  and  preconditioned  Krylov  iterative  linear  solvers.  We  then  make  recommen¬ 
dations  for  robust  and  efficient  approaches  to  obtain  steady-state  solutions  for  RE 
under  a  range  of  conditions. 


Notation 

Roman  Letters 

A  accumulation  term  for  pressure  head  form  of  RE 

A  accumulation  term  contribution  to  Jacobian 

J  Jacobian 

C  scaling  factor  for  Dirichlet  boundary  conditions 

F  nonlinear  function  for  DAE  formulation,  semi-discrete 

G  nonlinear  function  for  DAE  formulation,  discrete 

Ks  saturated  hydraulic  conductivity 

surface  saturated  hydraulic  conductivity 
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Ke  effective  hydraulic  conductivity 

Nmax  maximum  number  of  iterations  for  nonlinear  solution  methods 

O d  discrete  spatial  operator 

R  discrete  spatial  operator  and  source  term 

Se  effective  saturation 

T  extent  of  temporal  domain 

X d  extent  of  spatial  domain  along  Xd  axis 

f  source  term  for  the  aqueous  phase  evaluated  at  cell  centers 

/  source  term  for  the  aqueous  phase 

g  gravitational  vector 

g„  unit  vector,  g/||g|| 

kt  intrinsic  permeability 

kr  relative  permeability 

mv  parameter  for  VGM 

n  unit  outward  normal  for  D 

ne  total  number  of  nodes 

nv  parameter  for  VGM 

nxd  number  of  nodes  along  Xd  axis 

p  pressure  of  the  aqueous  phase 

t  time  coordinate 

u  mass  flux 

ub  Neumann  boundary  value 

ur  precipitation  rate 

y  variable  for  DAE  formulation 

y'  temporal  derivative  for  DAE  formulation 

Umin  lower  bound  for  solution  in  test  for  evaluation  error 

Umax  upper  bound  for  solution  in  test  for  evaluation  error 
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Greek  Letters 


At  time  step  for  'I' tc  methods 

A tmax  maximum  time  step  for  'htc  methods 

Axd  spatial  increment  in  Xa  direction 

Ay  Newton  increment 

T  boundary  of  physical  domain 

T e>  portion  of  T  for  which  Dirichlet  boundary  conditions  are  set 

T tv  portion  of  T  for  which  Neumann  boundary  conditions  are  set 

physical  domain 
av  parameter  for  VGM 

f3  compressibility  of  the  aqueous  phase 

ec  switching  tolerance  in  hybrid  Newton-Picard  method 

es  sufficient  decrease  parameter  for  Armijo  line  search 

6  volume  fraction  of  the  aqueous  phase 

6r  residual  volumetric  water  content 

6S  saturated  volumetric  water  content 

A  line-search  scaling  factor 

A+  intermediate  scaling  factor  for  quadratic  line  search 

/x  viscosity  of  the  aqueous  phase 

g  density  of  the  aqueous  phase 

p  normalized  density  of  the  aqueous  phase 

(Jrnvn  bound  parameter  for  quadratic  line  search 

crmax  bound  parameter  for  quadratic  line  search 

t  TTE  parameter  for  truncation  error  estimate  bound 

t/}  pressure  head 

,ipb  Dirichlet  boundary  value 
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initial  condition  for  i/j 


ip° 

Subscripts  and  Superscripts 

d  coordinate  axis  identifier  (subscript) 

i  cell  qualifier  along  x\  axis  (subscript) 

j  cell  qualifier  along  X2  axis  (subscript) 

l  global  identifier  for  solution  unknowns  (subscript) 

n  time  level  identifier  (superscript) 

sn  nonlinear  iteration  index  (superscript) 


Abbreviations 

Clock  wall  clock  time  (seconds) 

DAE  differential  algebraic  equation 

Feval  function  evaluation 

HAS  two-level  hybrid  additive  Schwarz  domain  decomposition  precon¬ 
ditioner 

Jeval  Jacobian  evaluation 

LF  linear  iteration  failure 

LI  linear  iteration 

LU-B  banded  LU  decomposition  from  LAPACK 

LU-S  sparse  LU  decomposition  from  SPOOLES  (version  2.2) 

NILS  Newton’s  method  with  a  quadratic  line  search 

NLF  nonlinear  iteration  failure 

NLI  nonlinear  iteration 

ODE  ordinary  differential  equation 
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PIH  Newton-Picard  hybrid  approach 

RE  Richards’  equation 

SER  switched  evolution  relaxation  \l/tc  approach 

Steps  attempted  iterations  for  method 

TTE  temporal  truncation  errror  d'tc  approach 

VGM  combined  van  Gennchten  and  Mnalern  p-S-k  relation 

p-S-k  pressure-saturation-relative  permeability  constitutive  relation 

'I'tc  pseudo-transient  continuation 


1  Introduction 

Variably  saturated  groundwater  flow  is  commonly  modeled  using  Richards’ 
equation  (RE)  along  with  a  set  of  constitutive  relations  describing  the  interde¬ 
pendence  among  fluid  pressures,  saturations,  and  relative  permeabilities  ( p-S - 
k  relations).  While  analytical  and  semi-analytical  approximations  for  variably 
saturated  flow  exist,  these  are  valid  for  limited  sets  of  auxiliary  conditions  and 
domains  [30].  As  a  result,  significant  effort  has  focused  on  developing  robust 
techniques  for  solving  RE  numerically  [9,  38,  29,  32,  27,  34,  36,  37].  Obtaining 
solutions  for  RE  for  many  realistic  physical  conditions  remains  a  challenge. 
Infiltration  problems  are  often  characterized  by  sharp  fronts  in  both  space 
and  time.  Steady-state  solution  is  often  nontrivial  as  well,  since  the  volume 
fraction  can  vary  steeply  for  problems  with  realistic  boundary  conditions  and 
heterogeneous  porous  media. 

It  is  the  nature  of  the  nonlinearities  introduced  through  standard  p-S-k  relations 
like  the  van  Genuchten  [35]  and  Mualem  [28]  models  that  accounts  for  the  ma¬ 
jority  of  the  difficulties  associated  with  solving  RE.  A  number  of  issues  associ- 
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ated  with  resolving  this  nonlinear  behavior  have  received  attention.  These  in¬ 
clude  solution  variable  transformations  [4,  36],  the  evaluation  oip-S-k  relations 
[33],  approximation  of  interface  conductivities  in  standard  low-order  spatial 
discretizations  [38,  27],  choice  of  dependent  variable  [37],  and  time  discretiza¬ 
tion  approach  [9,  31,  33,  19].  Since  the  majority  of  time  discretizations  are 
implicit,  both  transient  and  steady-state  problems  typically  lead  to  the  solu¬ 
tion  of  a  system  of  discrete  nonlinear  equations.  These  systems  can  be  quite 
large  and  difficult  to  solve,  especially  for  problems  with  heterogeneous  domains 
in  two  and  three  spatial  dimensions. 

The  nonlinear  system  solution  method  used  can  then  have  a  strong  impact  on 
the  overall  success  of  a  simulator  for  RE.  The  most  common  approaches  have 
been  Picard  iteration  [9,  10,  19]  or  a  variant  of  Newton’s  method  [34,  37,  21], 
Each  of  these  methods  has  its  strengths  and  weaknesses,  and  several  works 
have  compared  the  robustness  and  efficiency  of  Picard  and  Newton  approaches 
for  a  variety  of  problems  [29,  25,  27,  8].  Newton’s  method  combined  with  glob¬ 
alization  techniques  like  a  line  search,  reduction  of  time  step  (backtracking), 
or  the  use  of  Picard  iteration  to  obtain  an  initial  guess  has  proven  more  reli¬ 
able  and  efficient  than  Picard  iteration  for  several  problems  [29,  27].  However, 
both  Newton  and  Picard  approaches  have  been  shown  to  perform  poorly  un¬ 
der  certain  conditions.  The  solution  of  nonlinear  systems  can  be  particularly 
difficult  for  steady-state  problems  due  to  the  increased  difficulty  of  obtaining 
a  good  initial  guess  for  the  steady-state  solution.  In  addition,  the  scaling  of 
the  linearized  systems  is  typically  worse  [15],  and  there  is  no  longer  recourse 
to  reducing  time  steps  when  convergence  breaks  down  [29]. 

Pseudo-transient  continuation  (Ttc)  has  been  used  in  computational  fluid  dy¬ 
namics  for  some  time  to  obtain  steady-state  solutions  for  problems  in  which 
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Newton’s  method  with  standard  line-search  strategies  fails  [23,  24,  11],  Ttc 
combines  features  of  backward  Euler  time  integration  and  Newton’s  method 
with  the  idea  of  using  information  about  the  transient  physical  problem  to 
guide  selection  of  intermediate  iterates  to  the  steady-state  solution  [15].  Since 
it  uses  information  from  a  time-dependent  problem,  Ttc  can  be  more  robust 
than  standard  line-search  strategies,  while  incurring  less  computational  ex¬ 
pense  than  full  integration  of  the  transient  problem  [23]. 

In  this  work,  we  seek  effective  techniques  for  obtaining  steady-state  solutions 
to  RE.  We  examine  the  use  of  Ttc  methods  as  well  as  Newton’s  method 
combined  with  standard  globalization  techniques  for  steady-state  problems  in 
homogeneous  and  heterogeneous  domains.  We  investigate  the  methods’  per¬ 
formance  with  both  direct  and  preconditioned  Krylov  iterative  linear  solvers. 
We  then  make  recommendations  for  robust  and  efficient  approaches  to  obtain 
steady-state  solutions  for  RE  under  various  conditions. 


2  Background 

Historically,  Picard  iteration  has  been  the  most  common  nonlinear  solution 
method  used  for  RE  [9,  31].  It’s  appeal  can  be  traced  to  the  fact  that  it  is 
simple  to  implement,  since  it  does  not  involve  derivatives  of  terms  involving 
the  p-S-k  relations  and  produces  symmetric  linear  systems  [29].  While  Picard 
iteration  is  still  used  [10,  19],  Newton  and  inexact  Newton  approaches  have 
become  significantly  more  common  [14,  18,  21],  particularly  for  more  realistic 
multidimensional  problems.  An  inexact  Newton  approach  can  be  thought  of 
as  Newton’s  method  where  the  equation  for  the  Newton  update  is  satisfied 
only  approximately  [34],  This  may  be  due  to  the  fact  that  the  Newton  update 


is  solved  only  approximately  using,  for  example,  a  Krylov-type  iterative  linear 
solver  [18].  Or,  only  an  approximate  Jacobian  may  be  used  in  the  Newton 
update.  Common  approaches  for  approximating  the  full  Jacobian  are  to  use 
finite  differences  (a  numerical  Jacobian)  [14]  or  to  lag  its  evaluation  over  a 
number  of  iterates  (chord  or  modified  Newton’s  method)  [33]. 

While  these  approaches  may  be  conceptually  straightforward,  a  number  of 
subtleties  often  arise  [22],  For  instance,  the  choice  of  stopping  criteria  and  the 
associated  tolerances  for  both  the  nonlinear  and  linear  solvers,  the  frequency 
of  Jacobian  evaluations,  and  the  increment  used  in  evaluating  numerical  Ja- 
cobians  can  all  play  a  significant  role  in  the  success  of  an  approach  for  a 
given  problem  [29,  25,  27,  8,  34],  Since  the  convergence  of  Newton’s  method 
can  be  quite  sensitive  to  the  quality  of  initial  guess,  supplementing  a  New¬ 
ton  approach  with  strategies  like  a  line  search  and  backtracking  has  improved 
nonlinear  solver  performance  for  some  problems  [27].  Line-search  techniques 
typically  select  a  fraction  of  the  full  Newton  correction  to  update  the  current 
iterate  with  the  step  chosen  to  insure  that  the  updated  solution  represents  a 
sufficient  decrease  in  the  nonlinear  residual  [18,  11].  For  transient  problems, 
methods  commonly  adjust  the  time  step  based  on  nonlinear  solver  perfor¬ 
mance.  Specifically,  in  the  case  of  poor  nonlinear  solver  performance,  the  time 
step  is  often  reduced  [31,  33,  8].  This  can  improve  the  scaling  of  the  Jacobian, 
since  the  reciprocal  of  the  time  step  typically  appears  on  the  diagonal  [15], 
and  can  improve  the  quality  of  the  initial  guess  obtained  from  the  solution 
at  the  previous  time  step  [29].  Another  strategy  designed  to  address  Newton 
approaches’  sensitivity  to  initial  guess  is  a  hybrid  Newton-Picard  algorithm 
which  uses  a  Picard  iteration  initially  until  a  certain  level  of  convergence  is 
reached  and  then  switches  to  Newton’s  method  [29,  8]. 
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Like  the  hybrid  Newton-Picard  algorithm,  Ttc  attempts  to  combine  the  rapid 
convergence  of  Newton’s  method  near  the  steady-state  solution  with  a  more 
robust  iteration  far  from  the  solution.  Ttc  employs  a  time-stepping  approach 
when  the  iteration  is  far  from  the  steady  state,  and  thus  exploits  the  rela¬ 
tionship  between  the  nonlinear  system  for  the  steady-state  problem  and  the 
initial  value  problem  from  which  it  is  derived.  The  time-stepping  algorithm 
used  in  Ttc  has  properties  similar  to  forward  Euler  with  simple  heuristics, 
which  makes  it  both  unstable  and  inaccurate  as  a  time-integration  method  for 
stiff  ordinary  differential  equations  (ODE’s).  Nevertheless,  Ttc  can  effectively 
maintain  certain  important  qualities  of  the  transient  solution  of  some  prob¬ 
lems  such  as  enforcing  a  CFL  condition,  and  has  been  shown  to  be  effective  in 
reaching  the  steady  state  for  models  with  discontinuous  transient  phenomena 
such  as  shocks.  While  RE  does  not  exhibit  shocks  under  physically  realistic 
choices  of  parameters,  it  does  exhibit  sharp  fronts  due  to  the  nonlinearities  as 
well  as  steep  gradients  due  to  discontinuities  in  the  porous  media  properties. 


3  Approach 


3. 1  Formulation 


RE  can  be  formulated  in  a  number  of  ways.  We  begin  with  an  expression  for 
the  conservation  of  mass  for  the  aqueous  phase  in  an  air-water  system  where 
the  solid  phase  is  assumed  immobile,  and  interphase  mass  transfer  is  neglected. 
Combining  this  with  the  standard  extension  of  Darcy’s  law  to  variable  satu- 
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rated  conditions  [26]  gives 


djpO) 

dt 


v  •  pKe{yi>  -  pgu)  +  f 


in  ff  x  [0,  T] 


(1) 


with 

e  =  £>0e/3(p-Po) 
P  —  q/ Qo 


A 

_  goJM^i 

Ke  =  kM)Ks 

Here  p,  g,  and  p  are  the  pressure,  density,  and  viscosity  of  the  aqueous  phase, 
6  is  the  volume  fraction,  and  /  is  a  source  term  for  the  aqueous  phase.  £>o  is  the 
density  at  p0,  (3  is  the  compressibility  of  the  aqueous  phase,  ^  is  the  pressure 
head,  and  g  is  a  vector  accounting  for  the  acceleration  of  gravity.  Ke  is  the 
effective  hydraulic  conductivity,  Ks  is  the  saturated  hydraulic  conductivity, 
and  ki  is  the  intrinsic  permeability  of  the  porous  medium,  C  1R2  is  the 
spatial  domain,  and  [0,  T)  is  the  temporal  domain. 

The  auxiliary  conditions  for  eqn  (1)  are  given  by 


^  =  on  Yd,  t  €  [0,  T]  (3) 

u  •  n  =  ub  on  Tjv,  t  G  [0,  T]  (4) 

■0  =  il>°  in  £2,  t  —  0  (5) 

where 

u  =  -p/ie(V0  -  pgu)  (6) 
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is  the  mass  flux,  06,  and  ub  are  boundary  condition  functions.  0°  is  the  initial 
condition,  and  n  is  the  unit  outward  normal  for  Q.  We  have  also  set  T  =  dVt  = 
r £>  u  r at  with  r £>  n  r at  =  0. 

The  steady-state  form  of  eqn  (1)  is 

-V  •  p/ie(V0  -  pgu )  =  /  in  Q  (7) 


with  auxiliary  data  given  by  eqns  (3)  and  (4)  without  temporal  dependence. 

We  use  non-hysteretic  forms  of  the  p-S-k  relations  of  van  Genuchten  [35]  and 
Mualem  [28]  (VGM).  For  0  <  0,  these  are  given  by 


Se 


(0  -  Or) 

(0S  ~  0r ) 

[i+Ki0irrm’' 

i  -  (i  -  sl/m*y 


1  2 


(8) 

(9) 

(10) 


where  Se  is  the  effective  saturation,  0r  is  the  residual  volumetric  water  content, 
0S  is  the  saturated  volumetric  water  content,  av  is  a  parameter  related  to  the 
mean  pore  size,  nv  is  a  parameter  related  to  the  uniformity  of  the  soil  pore- 
size  distribution,  and  rnv  =  1  —  1  jnv.  For  0  >  0,  the  porous  medium  is  fully 
saturated,  and  eqns  (8)-(10)  revert  to 


Se  =  l  (11) 

kr  =  1  (12) 


The  pressure  head  form  of  RE  is  obtained  from  eqn  (1)  by  applying  the  chain 
rule  to  the  left  hand  side 


dip 

~dt 


V  •  p/ie(V0  -  pgu)  +  f 
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3.2  Spatial  discretization 


We  consider  a  discretization  of  Q  =  [0,Xl]  x  [0,X2]  C  IR2  into  a  regular, 
orthogonal  grid  with  ne  =  nxl  ■  nx 2  nodes  with  Axd  =  Xd/{nxd  —  1),  for  d  = 
1,  2.  We  apply  a  cell-centered  finite  difference  approximation  to  the  right  hand 
side  of  eqn  (1)  and  write  for  a  cell  in  the  interior, 

d(Pi,j@i,j)  _  J, 
dt  id 

11/ .j  Odtij  fij 

i  =  2, . . . ,  nxl  -  1,  j  =  2, . . . ,  nx2  -  1 


(14) 

(15) 


For  a  cell  O,; j  in  the  interior,  the  discrete  spatial  approximation  is 


Od,i,j 


Aaq 


Ax2 


('ll’i+lj  '  i>i,j 

1  Pi+l/2,jdul 

TV-  /  A,j  -  A~lj  \ 

Pi-l/2,jJ^e,i-l/2,j  I  ^  Pi-l/2,j9u\  J 

/  V’ij+i  -  Aj 

Pi,j+l/2^e,i,j+l/2  ^ ^ - Pi,j+l/29u2 

v  ( Aj  ~  A, j-l  \ 

Pi,j-l/2^e,i,j-l/2  I  - ^7 - Pi,j-l/29u2  j 


(16) 


where  the  subscripts  in  gu  =  [gui,gu2]J  indicate  values  for  each  coordinate, 
and  quantities  with  a  1/2  subscript  denote  values  estimated  at  cell  interfaces. 


For  the  interface  values,  we  use  a  harmonic  average  for  saturated  hydraulic 
conductivity,  the  arithmetic  average  for  density,  and  upwind  the  relative  per¬ 
meability 


Pi+i/2,,-  =  [p(V> i+u)  +  1/2 


(17) 


(18) 


-^■e  ,i+ 1  /  2  ,j  ^r,i+ 1  /2 ,  j  K s  1/2  j 

K s,i+l/2,j  sfaj K Sji+lrf / “I-  -^s,2+l,j) 

I  if  >  Pi+l/2,j9ul 

h'r.i—  1  /2.j  \ 

kr('tpij)  otherwise 

The  corresponding  terms  along  the  other  coordinate  axis  are  defined  symmet¬ 
rically. 

The  physical  boundaries  are  located  at  the  nodes  (cell  centers),  so  specify¬ 
ing  Dirichlet  conditions  in  pressure  or  volume  fraction  is  straightforward.  For 
example,  at  a  cell  f \j  C  T0,  we  set 

Ctyij  ~  <,-)  =  0,  (19) 

where  C  is  a  scaling  factor  that  gives  the  boundary  equation  roughly  the  same 
scaling  as  the  interior  nodes. 

For  cells  along  Tjv,  we  use  linear  extrapolation  to  apply  the  flux  at  the  exterior 
(artificial)  cell  boundary  rather  than  the  cell  center.  This  approach  allows  us 
to  use  the  same  nodal  equation  at  the  boundary  node  as  we  do  at  the  interior 
points  for  C  If,  for  example,  flij  C  Tn,  we  set 

ul,l/2,j  =  ~  ul,3/2,j  (20) 

at  the  fictitious  left  cell  boundary.  Here  u  =  [mi,u2]T- 

3. 3  Temporal  approximation 

Since  we  consider  both  direct  solution  of  the  steady-state  problem  eqn  (7)  and 
integration  of  eqn  (1)  to  steady  state,  we  begin  by  applying  the  cell-centered 
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finite  difference  spatial  approximation  to  eqn  (13)  in  a  method  of  lines  (MOL) 


context 


A(^, 


hi  > 


dt 


Ri,j 


i  —  nxi  —  1,  j  =  2, . . . ,  nx2  -  1 


(21) 


with  the  Dirichlet  and  Neumann  boundary  conditions  given  by  eqns  (19)  and 

(20). 

The  semi-discrete  system  corresponding  to  eqn  (21)  can  be  written  as  a  set  of 
differential  algebraic  equations  (DAE’s) 

F(t,y,y')=0  (22) 


where  F  represents  a  set  of  equations  that  depend  on  time  t,  a  set  of  dependent 
variables  y,  and  a  set  of  first-order  derivatives  with  respect  to  time  of  these 
dependent  variables,  y' . 


A  variety  of  approaches  can  be  used  to  integrate  eqn  (22).  To  illustrate  the 
structure  of  the  nonlinear  system  for  transient  problems,  we  use  a  backward 
Euler  approximation  to  convert  eqn  (22)  to  a  fully  discrete  system  [13,  21]. 


G(tn+1,  yn+1 


1 

’  A tn+1 


(y"+1 


y”))  =  o 


(23) 


3.3.1  Nonlinear  solution  for  the  transient  problem 

At  each  time  level,  a  full  time  integration  approach  such  as  backward  Euler 
must  solve  the  nonlinear  system  eqn  (23).  A  general  Newton  iteration  for  eqn 
(23)  can  be  written 

jn+M„]  |Ays"+1}  =  -  |Gn+1’Sn|  (24) 
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where  {AySn+1}  =  {yn+1,Sn+1}  —  {yn+1,Sn}  and  sn  is  a  nonlinear  iteration 
index.  The  Jacobian,  J,  is  formed  by  differentiating  eqn  (23)  with  respect  to 
y.  For  eqn  (21),  we  can  write  J  as 


rn+l,s„ 


1 

^n+l,s„ 

+ 

'<9(Ay') n+1,Sn" 

'dOdn+1,Sn~ 

Atn+1 

dy 

.  dy 

(25) 


where  A  is  diagonal  with  [A]n  =  Atyij),  dAy'/dy  is  also  diagonal,  and 
dO d/ dy  will  be  banded  with  seven  non-zero  entries.  Here,  /  is  a  global  identifier 
corresponding  to  cell  f \j  with,  for  example,  /  =  (j  —  l)nxi+i  for  i  =  1, . . . ,  nx i, 
j  !)•••)  Tlx2  ■ 

For  unknowns  along  Y Dl  J  is  simply  (see  eqn  (19)) 


rn+l,sn 


J  LI 


=  C 


(26) 


3-4  J7  tc:  approximation 


Ttc  attempts  to  hnd  a  solution  to  eqn  (7)  by  integrating  eqn  (21)  to  steady 
state.  The  approach  is  straightforward  and  can  included  in  many  existing 
steady-state  or  transient  solvers  with  minor  modifications.  The  fully  discrete 
Ttc  system  can  be  obtained  from  eqn  (21)  by  first  applying  a  backward  Euler 
time  discretization  as  in  eqn  (23)  and  then  using  Newton’s  method  with  yn 
as  the  initial  guess. 


3-4-1  Solution  for  'Ytc  update 

While  applying  multiple  iterations  of  Newton’s  method  is  possible,  the  form 
of  Ttc  which  we  present  here  performs  only  a  single  Newton  update  for  eqn 
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(23)  [23].  The  resulting  equation  for  the  Ttc  iterate  is 


[J“]  {Ay“+1}  =  -  {R“}  =  fO,/iy"i}  +  {f”} 


(27) 


with  {Ayn+1}  =  {yn+1}  —  {yn}  since  only  one  Newton  iteration  is  performed 
and  {yn+1'°}  =  {yra}.  J  has  the  same  form  as  eqn  (25),  but  without  a  derivative 
of  the  accumulation  term 


;jn; 


Atn+ 1  ^ 


90/' 

dy 


(28) 


Note  that  when  At  is  small  enough  [Jra]  «  A  *+1  [An]  and  therefore 


Ay”+1  «  -A tn+l  [An]_1  {Rn} 


(29) 


which  is  the  update  corresponding  to  the  forward  Euler  method  applied  to  the 
ODE  form  of  RE. 


3.Jh2  'I He  step  selection 


Ttc  solves  a  series  of  problems  of  the  form  in  eqn  (27),  while  adapting  the  time 
step  Atn+1  based  on  the  intermediate  solution’s  behavior.  There  are  a  number 
of  common  strategies,  for  selecting  Atn+1,  including  the  switched  evolution 
relaxation  (SER)  [23,  15,  11] 


Atn+1 


A  tn 


Rnl 

||Rn|| 


(30) 


The  temporal  truncation  error  approach  (TTE)  attempts  to  control  the  time 
step  based  on  the  local  temporal  truncation  error  [23,  15,  11].  It  chooses  A tn+1 
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so  that 


(A tn+1)2  d2yi 

2(TT 1# 


<  r 


(31) 


for  each  component  of  the  solution  yi  and  some  constant  r.  Here,  we  estimate 
d2yi/dt 2  at  tn  by  [11] 


d  2yi 
dt2 


(tn) 


2 

A  tn  +  At”-1 


yf  -  Vi 
At” 


n— 1 


tT*  -  t/zn“2 


At' 


n— 1 


(32) 


For  both  SER  and  TTE,  we  also  enforce  an  upper  bound  A tmax  on  the  chosen 
time  step. 


3.5  Newton’s  method  for  the  steady-state  problem 


Following  eqn  (21),  we  can  write  a  cell-centered  finite  difference  spatial  ap¬ 
proximation  of  eqn  (7)  to  solve  the  steady-state  problem  directly 


Rij  =  0  (33) 

i  =  2,...,nxi  -  1,  j  =  2, . . . ,  nx2  -  1 

with  Dirichlet  and  Neumann  boundary  conditions  given  by  eqns  (19)  and  (20) 
without  temporal  dependence. 

The  Newton  update  for  eqn  (33)  is 


[jsn]  |Ays"+1} 

[J] 


{Rs’1} 

dOd 

dy 


(34) 
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3.5.1  Globalization  techniques 


A  number  of  globalization  strategies  are  commonly  used  to  address  the  sensi¬ 
tivity  of  Newton’s  method  to  initial  guess.  The  Armijo  line-search  technique 
scales  the  original  Newton  update  by  a  factor  chosen  to  take  a  step  as  close  to 
the  original  update  as  possible  while  insuring  a  sufficient  decrease  in  the  non¬ 
linear  residual  [18].  At  iteration  level  sn  for  eqn  (34),  the  Armijo  line  search 
can  be  written 

(1)  Solve  eqn  (34)  for  AySn+1,  and  set  A  =  A+  =  1 

(2)  While  ||R(ySn  +  AAy^+1)||  >  (1  -  esA)||R(ys")|| 

Choose  A+ 

if  A+  &min A ,  A-)-  @min A 

else  if  A_|_  >  crmax A,  A_j_  &max A 
A  =  A+ 

(3)  Set  {ySn+1}  =  {ys"}  +  A  {Ays™+1} 

where  es  is  a  parameter  controlling  the  amount  of  decrease  required  in  the 
nonlinear  residual  and  A  is  the  final  scaling  factor.  Here,  we  chose  A+  so  that 
it  minimized  a  three-point  parabolic  approximation  of  ||i?(ySn  +  AAyS7l+1)||  = 
/(A).  Bounds  on  A+  are  dictated  by  amrn  and  crmax-  The  details  of  this  approach 
can  be  found  in  Kelley  [22],  For  this  work  we  set  crmin  =  0.1  ,amax  =  0.55,  and 
es  =  10”4.  The  line  search  can  fail  if  A  becomes  too  small.  In  this  case,  the 
nonlinear  solver  is  said  to  have  failed  due  to  line-search  stagnation  [11], 

To  avoid  evaluation  errors  in  the  constitutive  relations  and  to  increase  the 
robustness  of  the  iterations,  we  also  include  a  test  to  make  sure  that  the 
solution  is  within  broad,  physically  relevant  bounds. 
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(1)  Solve  eqn  (34)  for  AySn+1,  and  set  A  =  1 

(2)  While  y/n  +  AA y/n+1  £  [ymin,  ymax]  V  l 

A  =  A/2 

(3)  Set  {ySn+1}  =  {ys™}  +  A  {Ays™+1} 

In  the  numerical  experiments  presented  below,  we  use  an  interval  ymin  = 
-100  [m],  ymax  =  100  [m\. 

3.6  Picard  iteration 

Picard  iteration  has  been  used  widely  in  the  solution  of  RE  [10,  19].  In  brief, 
a  Picard  linearization  of  eqn  (33)  can  be  formulated  as  the  Newton  update  in 
eqn  (34)  but  with  an  approximate  Jacobian  in  which  the  coefficient  derivatives 
dkr/dip  and  dp/ dip  are  omitted  from  dOd/dy  [29,  25].  The  performance  of 
Picard  iteration  and  Newton  approaches  have  been  compared  in  several  works 
[29,  25,  27], 

In  many  cases,  Newton’s  method  with  line  search  has  proven  more  robust 
than  a  straightforward  application  of  Picard  iteration  [27].  However,  a  hybrid 
Newton-Picard  algorithm  has  also  been  suggested  for  reducing  the  sensitivity 
of  Newton’s  method  to  the  quality  of  the  initial  guess  [29,  8].  This  approach 
performs  an  initial  number  of  iterations  for  eqn  (34)  using  the  Picard  approx¬ 
imation  for  J  and  then  switches  to  a  Newton  update  with  the  full  Jacobian 
[29,  8].  The  motivation  for  this  approach  is  that  the  Picard  iterations  are,  in 
general,  cheaper  than  their  Newton  counterparts. 

Ideally,  the  switch  to  Newton’s  method  is  chosen  so  that  (1)  a  minimal  number 
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of  Picard  iterations  are  performed;  and,  (2)  upon  switching,  ySn  is  sufficiently 
close  to  the  true  to  solution  that  the  asymptotic  convergence  rate  for  the 
Newton  updates  is  realized.  There  are  several  possible  criteria  for  determining 
the  crossover  from  Picard  to  Newton  iteration,  including  ||Ays"||  <  ec  [8]. 
Alternatively,  one  can  base  the  switch  on  sufficient  decrease  in  the  initial 
residual, 

||RSn+1||  <  ec||R°||  (35) 

We  use  eqn  (35)  in  the  numerical  results  presented  below. 

3.7  Linear  system  solution 

We  test  five  algorithms  for  solving  the  linear  systems  arising  in  the  nonlinear 
iteration.  As  a  direct  solver,  we  use  both  the  banded  LU  decomposition  from 
LAPACK  (LU-B)  [1],  as  well  as  the  implementation  of  LU  from  the  SPOOLES 
package  (LU-S)  [2,  3]  in  PETSc  [6,  7,  5].  We  also  test  the  iterative  method 
BiCGstab  using  ILLI  preconditioning  with  zero  fill  from  PETSc  (BiCGstab- 
ILU),  and  a  two-level  hybrid  additive  Schwarz  domain  decomposition  method 
(BiCGstab- HAS)  [17]. 

4  Results 

In  the  following  sections  we  will  present  results  of  several  numerical  experi¬ 
ments  and  compare  the  behavior  of  Ttc  with  Newton’s  method  and  a  hybrid 
Newton-Picard  iteration.  First  we  summarize  the  test  problems  on  which  the 
numerical  experiments  were  carried  out. 
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4-1  Test  problems 


4-1.1  Problem  1:  infiltration  example 

The  first  test  problem  is  a  relatively  simple  one- dimensional  example  and  pro¬ 
vides  a  case  where  each  of  the  methods  should  perform  well.  We  consider  each 
test  case  as  a  transient  problem  with  a  corresponding  steady-state  solution. 
The  first  example  simulates  infiltration  in  vertical  domain  =  [0,  Xfi\.  Initial 
conditions  for  the  infiltration  are  set  to  static  equilibrium  with  the  water  table, 
located  at  the  bottom  of  the  domain.  Constant  Dirichlet  boundary  conditions 
are  set  at  the  top  so  that  the  steady-state  solution  contains  both  saturated 
and  unsaturated  regions.  Table  1  summarizes  the  relevant  physical  parameters 
and  auxiliary  conditions. 

4-1.2  Problem  2:  hillslope  example 

The  spatial  domain  for  the  second  problem,  =  [0,  x  [0,  X2],  is  illustrated 
in  Figure  1.  The  temporal  domain  is  t  G  [0,  T]  with  T  chosen  so  that  the 
solution  is  at  steady  state.  The  relevant  physical  parameters  and  auxiliary 
conditions  for  Problem  2  are  given  in  Tables  2-4.  The  log  of  the  saturated 
hydraulic  conductivity  is  given  in  Figure  2.  Note  that  the  domain  is  rotated 
45  degrees  to  simulate  a  simple  hillslope.  The  domain  has  block-heterogeneous 
medium  properties  ranging  from  clay  to  sand. 

The  boundary  and  initial  conditions  are  configured  to  reflect  an  imperme¬ 
able  bedrock  at  the  Xfi  boundary,  no  flow  at  the  A"]4-  boundary,  and  static, 
saturated  equilibrium  along  Xf .  The  surface  boundary  condition  at  Xfi  is 
a  nonlinear  flux  boundary  condition  that  simulates  infiltrating  precipitation 
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Table  1 


Fluid  and  domain  properties  for  Problem 

Variable  Value  Units 

9ui 

-1.0 

[-] 

||g|| 

7.321xl010 

[m/d2] 

go 

998.2 

CO 

.■Jg, 

P 

6.564xlO-20 

[m  ■  d/kg] 

Po 

0 

[kg /m  ■  d] 

Xx 

10 

[m] 

i/jb(xi  =  0) 

0 

[m] 

ij)b(x  i  =  10) 

-0.05 

[m] 

ip°(xi) 

-xi 

[m] 

nv 

4.264 

[-] 

QLy 

5.470 

[m-1] 

Ks 

5.040 

[m/d] 

until  the  surface  becomes  saturated.  After  saturated  conditions  are  reached  at 
the  surface,  the  outward  normal  flux  increases  to  zero  and  becomes  positive  as 
the  pressure  rises  above  atmospheric  pressure.  This  condition  reflects  a  well- 
drained  surface  that  permits  very  little  ponding  at  the  surface.  The  boundary 
conditions  can  be  summarized  as  follows. 


iix-=x-1 

-o 


(36) 

(37) 

(38) 


ur,  i/j  <  0 
ur  +  Kg'i/j,  i/j  >  0 


(39) 


where  ur  is  the  precipitation  rate  and  Kss  is  the  saturated  conductivity  at  the 
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surface.  For  the  numerical  experiments  presented,  Kss  =  10  \m/d] 

The  steady-state  solution  of  the  volume  fraction  is  given  in  Figure  3.  This  so¬ 
lution  was  computed  by  integrating  the  transient  problem  to  T  =  1700  [d]  at 
which  time  ||-R||oo  <  10~5.  The  solution  was  obtained  using  a  DAE  integrator 
with  relative  and  absolute  integration  tolerances  of  10-4  [20,  21],  The  tran¬ 
sient  solution  exhibits  steep  moving  fronts  throughout  the  domain  while  the 
equilibrium  solution  shown  maintains  a  number  of  steep  moisture  gradients 
due  to  the  layering  of  the  media. 

Fig.  1.  Problem  domain  D 


X2 


X+  (X1,X2) 


xr 


(0,0)  x2 


X 


+ 

1 


*-  X\ 


4-2  Numerical  experiments 

For  both  test  problems,  we  performed  several  numerical  experiments  on  a 
series  of  grids.  On  each  grid,  we  ran  Newton’s  method  with  a  quadratic  line 
search  (NILS),  the  Newton-Picard  hybrid  approach  (PIH),  as  well  as  two  Ttc 
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Table  2 


Media  distribution  for  Problem  2 


Medium  Type 

PI 

P2 

PI 

PI 

PI 

PI 

PI 

P2 

P2 

P2 

Left  [to] 

0 

0 

0.7 

1.4 

2.9 

7.1 

8.6 

2.9 

6.4 

9.3 

Bottom  [m] 

0 

0.1 

0.1 

0.1 

0.2 

0.1 

0.3 

0.1 

0.2 

0.1 

Right  [to] 

10 

0.71 

1.4 

2.9 

6.4 

9.3 

9.3 

7.1 

7.1 

10 

Top  [to] 

0.1 

1 

0.8 

0.3 

0.3 

0.3 

0.8 

0.2 

0.3 

1 

Medium  Type 

P2 

P2 

P2 

P2 

P2 

P3 

P3 

P3 

P4 

PI 

Left  [to] 

3.6 

1.4 

0.7 

1.4 

7.1 

1.4 

2.9 

5.0 

1.4 

5.0 

Bottom  [to] 

0.4 

0.6 

0.8 

0.7 

0.7 

0.4 

0.7 

0.4 

0.3 

0.5 

Right  [to] 

5.0 

8.6 

9.3 

2.9 

8.6 

3.6 

7.1 

8.6 

8.6 

8.6 

Top  [to] 

0.6 

0.7 

1 

0.8 

0.8 

0.6 

0.8 

0.5 

0.4 

0.6 

Table  3 

Media  properties 

for  Problem  2 

Medium  Type 

0s 

9r 

nv 

OLy 

[to.-1] 

Ks 

[m/d] 

PI 

0.41 

0.07749 

2.090 

0.244 

1.10808x10“ 

5 

P2 

0.40 

0.03120 

4.264 

5.470 

5.04000x10° 

P3 

0.39 

0.03822 

2.370 

0.478 

1.80100x10“ 

3 

P4 

0.39 

0.02691 

3.264 

0.244 

4.04000x10° 

methods,  SER  and  TTE,  until  ||i2||2  <  (H-R0II2  +  1)10“5.  The  initial  guess  was 
taken  to  be  the  initial  conditions  for  the  corresponding  transient  problem.  For 
the  NILS  and  PIH  calculations,  we  enforced  a  maximum  number  of  nonlinear 
iterations  ,sn  <  Nrnax  =  1000,  while  the  maximum  number  of  steps  for  the 
SER  and  TTE  runs  was  n  <  Nmax  =  5000.  The  maximum  number  of  line 
searches  allowed  was  1000  and  the  sufficient  decrease  parameter  for  Newton 
line  search  was  es  =  10“4.  The  PIH  approach  used  a  value  of  ec  =  10“2  in 
its  switching  strategy,  eqn  (35).  The  Ttc  methods  used  an  initial  time  step 


25 


y  [m] 


Table  4 


Fluid  and  domain  properties  for  Problem  2 


Variable 

Value 

Units 

9ui 

-0.7071 

[-] 

9ui 

-0.7071 

[-] 

g 

7.321xl010 

[m/d2] 

QO 

998.2 

GO 

P 

6.564xlO-20 

[m  ■  d/kg] 

Po 

0 

\kg/m  ■  d] 

10 

[m] 

x2 

1 

[, m ] 

ur 

-0.4 

[m/d] 

Fig.  2.  log(AT)  for  Problem  2 
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Fig.  3.  Steady-state  volume  fraction  for  second  test  problem 


x  [m] 


of  At0  =  4.050  x  10  5,  a  maximum  time  step  A tmax  =  105,  and  the  TTE 
algorithm  used  r  =  1.0.  The  relevant  parameters  are  summarized  in  Table  5. 

The  linear  systems  for  Problem  1  were  solved  using  the  LU-B  decomposition 
from  LAPACK  with  each  method.  For  Problem  2,  we  ran  the  numerical  exper¬ 
iments  for  each  of  the  four  steady-state  solution  methods  with  the  LU-S  direct 
solver  as  well  the  two  iterative  methods  (BiCGstab-ILU  and  BiCGstab-HAS). 
A  relative  residual  test  was  used  for  BiCGstab  with  a  tolerance  of  10 ”6.  The 
maximum  number  of  linear  iterations  allowed  was  2000  except  where  noted. 
The  simulations  were  performed  on  a  Pentium  4  (2.53  Ghz)  workstation  with 
256  Mbytes  of  RAM  running  Redhat  Linux  7.3.  The  elapsed  time  for  the 
simulations  was  recorded  using  the  ANSI  C  clock  and  time  intrinsics. 

The  results  for  Problem  1  on  spatial  grids  of  size  ne  =41-161  are  presented  in 


27 


Table  5 


Numerical  methods  summary 


NILS 

PIH 

Nmax 

1000 

1000 

es 

10~4 

- 

ec 

- 

10-2 

TTE 

SER 

Nmax 

5000 

5000 

At0 

4.050xl0-5 

4.050xl0-5 

^tmax 

105 

105 

T 

1.0 

Table  6.  The  labels  are  as  follows: 

•  Jeval-Jacobian  evaluation 

•  Feval- function  evaluation 

•  Steps-attempted  iterations  for  method 

•  NLI-nonlinear  iteration 

•  NLF-nonlinear  iteration  failure 

•  Ll-linear  iteration 

•  LF-linear  iteration  failure 

•  Clock-wall  clock  time  (seconds) 

•  X-failed  to  converge 

The  simulations  for  Problem  1  were  small  and  involved  a  homogeneous  medium. 
As  a  result,  we  expected  the  methods  to  have  little  difficulty  in  its  solution. 
Both  the  NILS  method  and  the  PIH  iteration  performed  well,  even  though 
the  NILS  method  employed  a  number  of  line  searches  on  each  spatial  grid. 


Table  6 


LU-B  runs  for  Problem  1 


Jeval 

Feval 

Steps 

NLI 

NLF 

LI 

LF 

LS 

Clock 

Tle  =  41 

NILS 

30 

299 

30 

30 

0 

0 

0 

102 

0.05 

PIH 

14 

33 

14 

14 

0 

0 

0 

1 

0.01 

TTE 

51 

103 

51 

0 

0 

0 

0 

0 

0.06 

SER 

83 

167 

83 

0 

0 

0 

0 

0 

0.04 

ne  =  81 

NILS 

26 

253 

26 

26 

0 

0 

0 

88 

0.04 

PIH 

15 

35 

15 

15 

0 

0 

0 

1 

0.03 

TTE 

104 

209 

104 

0 

0 

0 

0 

0 

0.14 

SER 

168 

337 

168 

0 

0 

0 

0 

0 

0.09 

ne  =  161 


NILS 

25 

329 

25 

25 

0 

0 

0 

128 

0.08 

PIH 

16 

35 

16 

16 

0 

0 

0 

0 

0.06 

TTE 

207 

410 

202 

0 

0 

0 

0 

0 

1.43 

SER 

353 

705 

351 

0 

0 

0 

0 

0 

0.33 

PIH  took  roughly  half  as  many  steps  as  NILS  and  needed  only  one  line  search 
once  it  switched  to  the  Newton  iteration,  indicating  that  the  initial  Picard 
iterations  were  more  successful  in  advancing  the  intermediate  iterates  closer 
to  the  steady-state  solution  than  Newton’s  method  with  a  line  search  alone. 
The  Ttc  methods  also  converged  to  the  correct  solution,  but  required  signifi¬ 
cantly  more  steps  than  either  the  NILS  or  PIH  approaches.  To  illustrate  the 
methods’  performance,  Figure  4  shows  the  residual  history  for  NILS  and  PIH. 
Both  methods  achieved  quadratic  convergence  as  they  approached  the  root. 


29 


Problem  2  was  more  challenging  than  Problem  1  clue  to  its  dimensionality,  het¬ 
erogeneous  medium,  and  the  nonlinear  boundary  condition  at  the  surface.  Re¬ 
sults  of  the  numerical  experiments  with  LU-S  and  BiCGstab  are  presented  in 
Tables  7-9.  The  increased  wall  clock  times,  iteration  counts,  and  line  searches 
reflect  the  added  difficulty  of  Problem  2.  NILS  converged  for  every  linear  solver 
and  spatial  grid  in  Tables  7-9.  It  required  a  similar  number  of  nonlinear  itera¬ 
tions  and  line  searches  for  a  particular  grid,  regardless  of  the  linear  solver  used. 
Unlike  NILS,  PIH  failed  for  the  cases  considered.  With  one  exception  where 
it  failed  in  the  initial  linear  solve,  PIH  did  not  reduce  the  original  nonlinear 
residual  sufficiently  to  satisfy  eqn  (35)  and  switch  to  NILS.  In  these  cases,  it 
exhausted  the  allowed  number  of  nonlinear  iterations. 

Both  the  Ttc  methods  converged  for  each  spatial  grid  and  linear  solver  com¬ 
bination  in  Tables  7-9,  with  TTE  consistently  between  2  and  5  times  faster 
than  SER.  The  run  times  for  NILS  and  TTE  were  similar  for  most  of  the  sim¬ 
ulations.  NILS  was  more  efficient  with  the  LU-S  solver,  particularly  for  the 
81  x  81  grid.  On  the  other  hand,  TTE  was  more  efficient  for  the  BiCGstab 
calculations.  The  difference  in  run  times  increased  for  the  HAS  preconditioner, 
where  TTE  was  3  and  1.4  times  faster  than  NILS  on  the  41  x  41  and  81  x  81 
grids  respectively. 

The  results  also  indicate  a  difference  in  the  way  NILS  and  TTE  behaved  as 
the  spatial  grid  was  refined.  Namely,  NILS  required  roughly  the  same  number 
of  iterations  to  converge  regardless  of  spatial  grid  or  linear  solver,  while  the 
number  of  steps  taken  by  TTE  grew  noticeably  as  ne  increased.  To  investigate 
the  performance  of  TTE,  we  ran  an  additional  set  of  calculations  where  r  was 
set  to  1CU3  and  scaled  by  ne  for  each  grid.  The  corresponding  r  values  ranged 
from  0.121  on  the  11  x  11  grid  to  25.9  on  a  grid  with  ne  =  161  x  161.  As  a 
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result,  the  TTE  time  step  selection  was  slightly  more  conservative  on  coarser 
grids  than  the  original  choice  of  r  =  1  and  was  more  aggressive  on  the  finer 
grids.  Table  10  contains  the  results  for  TTE  with  LU-S  and  BiCGstab  for 
ne  =  11  x  11  to  ne  =  161  x  161  as  well  as  the  results  for  NILS  on  the  finest 
grid. 

With  the  use  of  a  scaled  r,  Steps  for  TTE  was  roughly  the  same  for  ne  =  11x11 
to  ne  =  81  x  81  and  each  linear  solver.  There  was  still,  however,  an  increase 
in  the  number  of  steps  for  the  finest  grid,  particularly  for  the  BiCGstab-ILLI 
solver.  TTE  with  the  scaled  r  took  slightly  longer  than  the  original  r  =  1 
calculations  on  the  coarser  grids,  but  reduced  the  total  simulation  time  as  the 
grids  were  refined.  There  was  a  corresponding  improvement  for  TTE  on  the 
larger  spatial  grids  with  all  three  solvers.  The  computational  effort  for  NILS 
and  TTE  with  r  scaled  was  essentially  the  same  for  LU-S  except  for  the  finest 
grid,  where  NILS  was  1.6  times  faster.  For  BiCGstab- ILLI  and  BiCGstab- 
HAS,  TTE  with  r  scaled  was  twice  as  fast  as  the  simulations  with  r  =  1  on 
the  ne  =  81  x  81  grid. 

We  also  note  that  BiCGstab-ILU  did  not  perform  well  with  the  ne  =  161  x 
161  system.  Simulations  with  each  of  the  steady-state  solution  methods  and 
BiCGstab-ILU  typically  required  significantly  less  computational  effort  than 
BiCGstab-HAS  for  the  coarser  grids.  However,  for  ne  =  161  x  161  the  total 
number  of  steps,  linear  iterations,  and  simulation  time  increased  significantly 
for  TTE  while  NILS  failed  after  both  2000  and  10000  linear  iterations. 

As  an  example  of  the  progress  of  the  Ttc  iterations  in  comparison  to  the 
Newton  and  Newton-Picard  iterations,  we  graphed  the  history  of  ||i?||2  versus 
iteration  for  SER  on  the  ne  =  81x81  grid  with  the  LLI-S  linear  solver  in  Figure 
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Fig.  4.  NILS/PIH  residual  history,  Problem  1 


5.  As  it  neared  the  steady-state  solution,  the  method’s  convergence  steepened. 
This  reflects  the  convergence  of  the  Newton  iteration  that  Ttc  reduces  to  near 
the  root.  Unlike  the  globalized  Newton  iteration,  however,  Ttc  does  not  enforce 
a  decreasing  sequence  of  residuals,  and  Figure  5  shows  increased  residuals 
at  some  points  in  the  SER  calculations.  The  iteration  history  of  TTE  was 
more  extreme,  often  oscillating  wildly  before  nearing  the  steady-state  solution. 
Figure  6  shows  the  TTE  calculation  with  r  =  1  on  the  ne  =  81  x  81  grid. 
For  some  problems,  it  has  been  found  that  performing  multiple  sub-iterations 
for  eqn  (27)  can  improve  convergence  [23,  15].  However,  we  investigated  this 
alternative  for  the  SER  and  TTE  updates  for  the  second  test  problem  and 
found  no  improvement  for  either  approach. 
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Table  7 


LU-S  runs  for  Problem  2 


Jeval 

Feval 

Steps 

NLI 

NLF 

LI 

LF 

LS 

Clock 

ne  =  11x11 

NILS 

83 

3901 

83 

83 

0 

83 

0 

1806 

1.5 

PIH 

1000 

6173 

1000 

1000 

1 

1000 

0 

0 

X 

TTE 

49 

99 

49 

0 

0 

49 

0 

0 

0.26 

SER 

145 

291 

145 

0 

0 

145 

0 

0 

1.65 

ne  =  21  x  21 

NILS 

132 

7057 

132 

132 

0 

132 

0 

3291 

2.55 

PIH 

1000 

2239 

1000 

1000 

1 

1000 

0 

0 

X 

TTE 

91 

183 

91 

0 

0 

91 

0 

0 

2.59 

SER 

379 

759 

379 

0 

0 

379 

0 

0 

7.5 

ne  =  41  x  41 

NILS 

106 

4725 

106 

106 

0 

106 

0 

2181 

9.92 

PIH 

1000 

2885 

1000 

1000 

1 

1000 

0 

0 

X 

TTE 

173 

346 

172 

0 

0 

173 

0 

0 

12.23 

SER 

866 

1732 

865 

0 

0 

866 

0 

0 

53.29 

ne  =  81  x  81 

NILS 

149 

7369 

149 

149 

0 

149 

0 

3412 

56.44 

PIH 

1000 

3449 

1000 

1000 

1 

1000 

0 

0 

X 

TTE 

450 

888 

437 

0 

0 

450 

0 

0 

164.43 

SER 

2068 

4135 

2066 

0 

0 

2068 

0 

0 

709.67 

5  Discussion 


We  began  with  the  goal  of  identifying  effective  methods  for  the  steady-state 
solution  of  RE.  To  this  end,  we  performed  several  numerical  experiments  for 
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Table  8 


BiCGstab-ILU  runs  for  Problem  2 


Jeval 

Feval 

Steps 

NLI 

NLF 

LI 

LF 

LS 

Clock 

ne  =  11x11 

NILS 

83 

3901 

83 

83 

0 

1297 

0 

1806 

1.58 

PIH 

1000 

6173 

1000 

1000 

1 

7371 

0 

0 

X 

TTE 

49 

99 

49 

0 

0 

210 

0 

0 

0.3 

SER 

145 

291 

145 

0 

0 

445 

0 

0 

1.77 

ne  =  21  x  21 

NILS 

134 

7205 

134 

134 

0 

6895 

0 

3360 

4.5 

PIH 

1000 

2239 

1000 

1000 

1 

19419 

0 

0 

X 

TTE 

98 

197 

98 

0 

0 

970 

0 

0 

2.08 

SER 

378 

757 

378 

0 

0 

2559 

0 

0 

7.14 

ne  =  41  x  41 

NILS 

111 

4863 

111 

111 

0 

11153 

0 

2239 

17.01 

PIH 

1000 

2893 

1000 

1000 

1 

38981 

0 

0 

X 

TTE 

177 

354 

176 

0 

0 

2215 

0 

0 

13.49 

SER 

865 

1730 

864 

0 

0 

11693 

0 

0 

54.44 

ne  =  81  x  81 

NILS 

145 

7127 

145 

145 

0 

31799 

0 

3301 

148.57 

PIH 

1 

3 

1 

1 

0 

2000 

1 

0 

X 

TTE 

424 

847 

422 

0 

0 

7934 

0 

0 

136.65 

SER 

2068 

4135 

2066 

0 

0 

50243 

0 

0 

604 

Newton’s  method  with  two  globalization  techniques  (NILS  and  PIH)  as  well 
as  two  versions  of  d'tc  (SER  and  TTE).  Various  observations  can  be  made 
based  on  the  results,  which  reflect  a  range  of  difficulty  for  the  one  and  two- 
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Table  9 


BiCGstab-HAS  runs  for  Problem  2 


Jeval 

Feval 

Steps 

NLI 

NLF 

LI 

LF 

LS 

Clock 

Tle  = 

11  x  11 

NILS 

84 

3959 

84 

84 

0 

4200 

0 

1833 

1.13 

PIH 

1000 

6173 

1000 

1000 

1 

24663 

0 

0 

X 

TTE 

49 

99 

49 

0 

0 

425 

0 

0 

0.38 

SER 

146 

293 

146 

0 

0 

1076 

0 

0 

1.99 

Tle  = 

21  x  21 

NILS 

131 

7059 

131 

131 

0 

46899 

0 

3293 

26.3 

PIH 

1000 

2239 

1000 

1000 

1 

51097 

0 

0 

X 

TTE 

88 

177 

88 

0 

0 

1628 

0 

0 

2.63 

SER 

379 

759 

379 

0 

0 

5799 

0 

0 

9.04 

Tie  = 

41  x  41 

NILS 

110 

4991 

110 

110 

0 

24471 

0 

2305 

69.02 

PIH 

1000 

3019 

1000 

1000 

1 

74767 

0 

0 

X 

TTE 

168 

336 

167 

0 

0 

4814 

0 

0 

23.93 

SER 

866 

1732 

865 

0 

0 

26241 

0 

0 

117.69 

Tie  = 

81  x  81 

NILS 

148 

7315 

148 

148 

0 

30762 

0 

3389 

395.15 

PIH 

1000 

3449 

1000 

1000 

1 

69639 

0 

0 

X 

TTE 

439 

867 

427 

0 

0 

11964 

0 

0 

288.42 

SER 

2068 

4135 

2066 

0 

0 

78706 

0 

0 

1377.61 

dimensional  variably  saturated  flow  problems  examined  in  this  work. 

The  use  of  a  line  search  or  hybrid  Picard  iteration  made  Newton’s  method 
more  robust  for  Problem  1.  PIH  required  roughly  half  the  iterations  of  NILS, 
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Table  10 


Runs  with  r  scaled  by  ne  for  Problem  2 


ne 

Jeval 

Feval 

Steps 

NLI 

NLF 

LI 

LF 

LS 

Clock 

LU-S 

TTE 

11  x  11 

114 

229 

114 

0 

0 

114 

0 

0 

1.57 

TTE 

21  x  21 

127 

255 

127 

0 

0 

127 

0 

0 

2.96 

TTE 

41  x  41 

125 

250 

124 

0 

0 

125 

0 

0 

8.16 

TTE 

81  x  81 

146 

281 

134 

0 

0 

146 

0 

0 

58 

TTE 

161  x  161 

248 

458 

209 

0 

0 

248 

0 

0 

667.07 

NILS 

161  x  161 

232 

12961 

232 

232 

0 

232 

0 

6043 

429.58 

BiCGstab-ILU 

TTE 

11  x  11 

113 

227 

113 

0 

0 

447 

0 

0 

0.69 

TTE 

21  x  21 

129 

259 

129 

0 

0 

1132 

0 

0 

2.36 

TTE 

41  x  41 

130 

260 

129 

0 

0 

1739 

0 

0 

9.72 

TTE 

81  x  81 

148 

289 

140 

0 

0 

3474 

0 

0 

65.85 

TTE 

161  x  161 

806 

1185 

378 

0 

0 

92232 

0 

0 

2330 

NILS 

161  x  161 

1 

3 

1 

1 

0 

10000 

1 

0 

X 

BiCGstab-HAS 

TTE 

11  x  11 

115 

231 

115 

0 

0 

979 

0 

0 

1.86 

TTE 

21  x  21 

125 

251 

125 

0 

0 

2173 

0 

0 

4.28 

TTE 

41  x  41 

131 

262 

130 

0 

0 

3782 

0 

0 

18.17 

TTE 

81  x  81 

152 

288 

135 

0 

0 

7192 

0 

0 

131.9 

TTE 

161  x  161 

283 

472 

188 

0 

0 

20806 

0 

0 

1386.07 

NILS 

161  x  161 

236 

13153 

236 

236 

0 

70484 

0 

6130 

3397 

36 


Fig.  5.  SER  residual  history,  Problem  2 


Fig.  6.  TTE  residual  history,  Problem  2 


step 


while  both  it  and  NILS  took  considerably  fewer  steps  than  either  SER  or 
TTE  to  reach  steady  state.  The  first  test  problem  was  intended  to  represent 
behavior  for  mild  problems  since  the  boundary  conditions  were  simple  and 


37 


the  conductivity  was  homogeneous.  The  resulting  linear  systems  were  also 
small  and  tridiagonal.  Under  these  ideal  conditions,  NILS  had  little  difficulty, 
and  PIH  was  able  to  speed  convergence  to  the  steady-state  solution.  The  Ttc 
methods  offered  no  benefit  in  this  situation. 

The  methods’  relative  performance  changed  though,  as  we  moved  to  more  re¬ 
alistic  conditions.  The  controlling  factor  in  their  performance  for  Problem  2 
was  the  difficulty  associated  with  solving  the  resulting  linear  systems  for  each 
method.  In  addition  to  nonlinear  boundary  conditions,  Problem  2  involved  a 
two-dimensional,  block- heterogeneous  domain  with  medium  properties  corre¬ 
sponding  to  sand  and  clay.  As  a  result,  the  linear  systems  for  Problem  2  were 
more  poorly  scaled  than  in  Problem  1.  On  the  first  step  of  the  81  x  81  grid  the 
condition  number  of  the  NILS  Jacobian  was  approximately  1019  and  after  150 
iterations  was  109  (condition  numbers  calculated  using  the  dgbtrf /dgbcon 
routines  from  LAPACK). 

Poor  scaling  of  the  Jacobians  manifested  itself  in  a  number  of  ways,  including 
the  added  work  required  by  BiCGstab  to  converge  with  both  preconditioners. 
In  general,  the  difficulty  of  solving  the  linear  systems  can  also  translate  into 
less  accurate  search  directions  as,  for  example,  one  is  forced  to  use  coarser 
linear  solver  tolerances  to  obtain  results  given  computational  limitations.  The 
performance  of  NILS  with  the  LLI-S  solver  illustrates  the  impact  of  the  Jaco¬ 
bian  scaling  on  direct  methods  and  the  result  of  inaccurate  search  directions. 
The  results  in  Table  7  indicate  that  NILS  was  the  most  efficient  approach  when 
LLI-S  was  used  for  the  LLI  decompositions.  However,  when  we  performed  the 
tests  with  the  LU  solver  from  PETSc  which  does  not  include  pivoting,  the 
solution  was  less  accurate,  and  NILS  either  exhausted  the  allowed  nonlinear 
iterations  or  stagnated  in  the  line  search. 
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The  PIH  approach  did  not  perform  as  it  was  intended  for  Problem  2  because 
the  initial  Picard  iterations  did  not  adequately  reduce  the  initial  residual. 
The  convergence  of  PIH  could  have  been  improved  by  relaxing  the  switching 
criterion  in  eqn  (35)  or  enforcing  a  maximum  number  of  Picard  iterations, 
so  that  it  would  have  switched  to  the  Newton  updates  before  exhausting  the 
allowed  iterations.  However,  this  would  be,  in  essence,  just  NILS  and  would 
have  masked  the  ineffectiveness  of  the  initial  Picard  iterations  for  Problem  2. 

In  contrast,  SER  was  robust  but  relatively  inefficient.  It  converged  for  each  of 
the  runs  and  linear  solvers,  but  was  significantly  slower  than  TTE  for  most 
cases.  The  two  most  efficient  methods  for  Problem  2  were  NILS  and  TTE. 
Their  relative  efficiency  was  dictated  largely  by  the  computational  expense 
associated  with  the  linear  solvers.  With  both  methods,  the  number  of  steps 
required  for  a  given  spatial  grid  was  similar  using  the  LLI-S  and  BiCGstab 
solvers.  However,  the  computational  effort  required  changed  significantly.  Us¬ 
ing  a  robust  direct  solver,  NILS  was  three  time  faster  than  TTE  with  r  =  1 
on  the  81  x  81  grid.  However,  it  required  four  times  as  many  total  linear  itera¬ 
tions  with  BiCGstab-ILU  and  three  times  as  many  using  BiCGstab-HAS.  As  a 
result,  it  was  9%  slower  using  BiCGstab-ILU  and  28%  slower  using  BiCGstab- 
HAS. 

The  total  number  of  iterations  for  NILS  was  fairly  consistent  on  the  spatial 
grids  considered.  On  the  other  hand,  TTE  with  r  =  1  did  not  scale  as  well 
when  the  spatial  grids  were  refined.  The  total  number  of  steps  required  to 
converge  increased  with  ne ,  making  TTE  more  expensive.  Varying  r  based  on 
ne  led  to  more  consistent  results  for  TTE  and  improved  its  performance  for 
finer  grids.  NILS  was  still  more  efficient  with  a  direct  solver  on  the  161  x  161 
grid.  On  the  other  hand,  TTE  was  between  2  and  3  times  faster  than  NILS  on 
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the  finer  grids  using  the  iterative  linear  solvers  and  converged  for  BiCGstab- 
ILU  on  the  161  x  161  grid  when  NILS  failed. 

The  value  of  a  number  of  parameters  like  ec,  es,  crmin,  crmax  and  the  linear  solver 
tolerances  for  BiCGstab  effected  each  of  the  methods  considered  to  some  de¬ 
gree.  The  performance  of  TTE  for  different  r  values  demonstrates  that  the 
the  impact  of  parameter  choices  can  be  significant  at  times.  Finding  a  suc¬ 
cessful  configuration  for  a  given  method  and  problem  usually  requires  some 
effort.  While  different  choices  of  parameter  values  may  lead  to  improved  per¬ 
formance,  the  basic  behavior  of  the  methods  was  the  same  for  the  range  of 
parameters  we  investigated.  NILS  was  preferable  for  the  simulations  with  a 
robust  direct  solver  like  LLI-S,  but  BiCGstab  performed  better  with  the  Ttc 
algorithms.  TTE  was  more  aggressive  than  SER  and  more  efficient  than  NILS 
when  BiCGstab  was  used  with  either  ILU  or  HAS  preconditioning. 

As  a  reference  point,  we  also  compared  the  performance  of  the  Ttc  methods 
to  time-accurate  integration  using  a  variable-order,  variable  step-size  DAE 
integrator  [20].  As  might  be  expected,  the  full  time  integration  was  the  most 
reliable  approach  for  obtaining  steady-state  solutions,  but  the  computational 
expense  incurred  was  from  5  to  25  times  greater  than  that  of  SER,  which  was 
similarly  robust. 

An  initial  guess  closer  to  the  final  solution  would  have  improved  the  perfor¬ 
mance  of  each  of  the  methods  considered.  Since  the  goal  of  PIH  and  the  Ttc 
methods  is  to  approximate  a  Newton  iteration  near  the  steady-state  solution, 
one  can  expect  the  advantage  of  NILS  to  increase  as  the  initial  guess  better 
approximates  the  steady-state  solution.  However,  we  used  static  equilibrium 
as  a  reasonable  initial  guess,  since  we  are  interested  in  evaluating  the  rneth- 
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ods  for  problems  where  good  initial  guesses  for  the  steady-state  solution  are 
difficult  to  obtain. 


The  range  of  spatial  grids  presented  for  the  numerical  experiments  was  rela¬ 
tively  coarse  and  the  resulting  linear  systems  were  small  to  moderate  in  size. 
Still,  the  Jacobians  were  ill-conditioned  for  Problem  2,  which  proved  to  be  a 
significant  test  for  the  direct  and  iterative  linear  solvers.  For  larger  systems 
arising  from  realistic  two  and  three-dimensional  problems,  we  can  only  ex¬ 
pect  the  difficulties  associated  with  the  linear  systems  to  become  more  severe. 
Moreover,  large  simulations  will  often  have  to  be  solved  in  parallel  to  obtain 
results  in  a  reasonable  timeframe.  For  the  numerical  experiments  presented 
here,  NILS  combined  with  the  LU-S  solver  was  the  most  efficient  approach  on 
each  spatial  grid.  However,  preconditioned  Krylov  methods  and  sparse  direct 
solvers  have  their  own  advantages  and  disadvantages  depending  on  the  prob¬ 
lem  and  architecture  [12,  16].  A  general  comparison  of  their  relative  merits  is 
beyond  the  scope  of  this  paper. 


6  Conclusions 

Our  numerical  experiments  for  Ttc  approaches  as  well  as  Newton’s  method 
with  various  globalization  techniques  lead  us  to  the  following  conclusions  and 
recommendations: 

•  For  problems  where  use  of  a  robust  direct  solver  is  feasible,  Newton’s  method 
with  a  line  search  is  the  most  efficient  approach  for  obtaining  steady-state 
solutions  to  RE. 

•  Using  an  initial  number  of  Picard  iterations  for  Newton’s  method  with  line 
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search  can  improve  performance  in  some  instances,  but  it  does  not  neces¬ 
sarily  lead  to  more  robust  performance  for  difficult  problems. 

•  Inexact  Newton  methods  with  standard  globalization  techniques  have  par¬ 
ticular  difficulty  when  the  Jacobian  is  poorly  scaled  due  to  factors  such  as 
heterogeneous  conductivity  fields. 

•  If  Newton’s  method  fails  or  performs  poorly  for  a  given  steady-state  prob¬ 
lem,  it  is  worth  examining  a  range  of  linear  solver  and  line-search  parameters 
before  abandoning  a  Newton  approach. 

•  Ttc  is  a  relatively  simple  approach  that  can  improve  the  efficiency  and 
robustness  of  existing  steady-state  solvers  for  RE  on  difficult  problems,  par¬ 
ticularly  if  iterative  linear  solvers  are  used. 
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